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ABSTRACT 

Toy Stars are gas masses where the compressibihty is treated without approximations 
but gravity is replaced by a force which, for any pair of masses, is along their line of 
centres and proportional to their separation. They provide an invaluable resource for 
testing the suitability of numerical codes for astrophysical gas dynamics. In this paper 
we derive the equations for both small amplitude oscillations and non linear solutions 
for rotating and pulsating Toy Stars in two dimensions, and show that the solutions 
can be reduced to a small number of ordinary differential equations. We compare 
the accurate solutions of these equations with SPH simulations. The two dimensional 
Toy Star solutions are found to provide an excellent benchmark for SPH algorithms, 
highlighting many of the strengths and also some weaknesses of the method. 
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1 INTRODUCTION 

A fundamental computational problem in astrophysics is the 
■ motion of a cloud of gas forming a protostar in an ambient 
medium which is typically lower in density by a factor of 
< 10^^. The ambient medium can therefore be treated as 
a vacuum to an excellent approximation. The simulation of 
. such a gas cloud is more difficult than many of the standard 
' battery of test problems considered in computational astro- 
physics. Exact, or very accurate, solutions for the motion of 
these gas clouds are rare and this creates further difficulties 
in assessing the accuracy of a numerical algorithm. 

In this paper we fill this gap in test cases by studying 
solutions and simulations for the class of models which we 
have called Toy Stars. They are, first and foremost, a tool 
for studying the accuracy of computational algorithms rele- 
vant to astrophysical gas dynamics but they also provide an 
interesting class of dynamical systems where exact solutions 
can be found. The fundamental aspects of Toy Stars and a 
variety of one dimensional solutions have been discussed by 
^onaghan & Price (2004). The key features are that they 
are gas masses where compressibility is included without ap- 
proximation, but gravity is replaced by a force derived from 
the potential 
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where v is a. constant, rrii is the mass of particle i and it is 
assumed that there are N masses. The equation of motion 
of particle j in the absence of other forces is 
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dt 

where M is the total mass and the origin is the cen- 
tre of mass. Remark ably, as first noted by Newton (see 
IChandrasekhaj ^9^, the particles move independently 
about the centre of mass of the particle system. Two parti- 
cles in a binary Toy Star system move on orbits given by the 
solutions of the following equation for the relative coordinate 
r = ri - r2 
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These solutions are closed Lissajous figures which include 
elliptical orbits. 

The equations of motion of a gaseous system in a Toy 
Star force are the acceleration equation 

dt p ^ ' 

where = Mv, P is the pressure and p the density, and 
the continuity equation 



dt ^ 



(5) 



If P = Kp^ this equation together with the continuity equa- 
tion is identical to the equation for small oscillations of wa- 
ter in a lake with paraboloidal bottom. This equation has 
been studied bv .Gold sbrough (,1930j and .Lamb (,193^'l. but 
the m ost important co ntributions were made bv lBalll jjlQel 
andlH^I^ ^^S) • In this paper we focus on the aspects 
of these equations which are most important for astrophys- 
ical problems, and generalise to the case where P = Kp"^ . 
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In the following we first consider the small amplitude 
oscillations of the Toy Star and compare the results to simu- 
lations using SPH. We then analyze the non linear motions 
for arbitrary 7 and show that the equations reduce to a 
small set of ordinary differential equations which can be in- 
tegrated with high accuracy. We then compare the results 
of these integrations with SPH simulations. 



2 STATIC STRUCTURE AND OSCILLATIONS 

If P = Kp'' the density of the static model is given by 
^2 ^ 
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where the radius re is 

" n2(7-l)' 
The mass M is given by 
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If the |v| < Cs, where Cs is the unperturbed speed of 
sound, the equations of motion can be linearized about the 
static structure. It is convenient to define 

R = p^-^ (9) 

and write R = R + rj where R is R calculated using the un- 
perturbed density p. The equations of motion then become 

'dt ^ " 
dr] 

We assume the time variation is e"'*, and write 
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If these expressions are substituted into the linearized equa- 
tions, V can be eliminated to get the following equation for 
D 
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Assuming separable solutions, they must be of the form 
C,{r) sinO or C,{t) cosO and the equation for is 

/d^C IdC s^C 
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The solutions of this equation determine the values of a. 
While this equation can be transformed to the equation for a 
Hypergeometric function it is more convenient to determine 
the solutions directly using expansions in series following the 
method of Frobenius. We thus take 



C(r) =X=^a„X", 
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where X — r /r^. It is convenient to replace o hy v according 
to 
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using the definition of re. If the series is substituted into the 
equation for we get the following recurrence relation for 
the coefficients 



afc+2 = flfe- 
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The indicial equation gives c — s and there is one so- 
lution with ao arbitrary and ai zero. Because the equation 
is second order there must be two arbitrary constants but 
because the solutions of the indicial equation differ by an 
integer (s is an integer) or are equal (s=0), the second ar- 
bitrary constant multiplies a solution containing In a;, and 
must be zero. The remaining series only converges if 
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where j is an integer and the associated value of a is denoted 
by Uj. The last term in the series for a given j is ajX-' . For 
numerical work we write the velocity in the form 



V = V cos {at) 
r) = D sin (at). 

From (Unj 
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Further details of the linear modes are given in Appendix 1X1 



V = 



-VD. 



3 SPH SIMULATIONS OF STATIC 

STRUCTURE AND LINEAR OSCILLATIONS 

Smoothed Particle Hydrodynamics (SPH) (for a recent re- 
view see Monaghan 2005i) is a Lagrangian particle method 
for solving the equations of fluid dynamics. Since a primary 
application of SPH is to self-gravitating gas in astrophys- 
ical systems, in most cases involving free boundaries. Toy 
Stars represent an ideal test of the algorithm's capabilities 
on these systems. Whil st the standard S PH algorithm (as de- 
scribed, for example in lMonaghaJl992ll has been well tested 
and benchmarked, we use the opportunity provided by the 
Toy Star solutions to benchmark more recent improvements 
to the algorithm. 

In particular we formulate the SPH equations from 
a variational principle such that the spatial variation of 
the smoothing length according to the density variation is 
accounted for self-consistently (Springcl & Hcrnauist 200j; 
lMonagharJl2002l: IPrice fc Monaghanll2004l: lPricdl2004D . We 
also use the Toy Stars to test a reversible time integration 
algorithm for SPH described in lMonagharJ (|200^. The spe- 
cific implementation of the SPH algorithm used for the test 
problems presented in this paper is described below. 



3.1 SPH implementation 

The SPH equations are formulated from a variational princi- 
ple which accounts for the spatial variation of the smoothing 
length with density. Prior to the force evaluation, the density 
is calculated via a direct summation over the particles which 
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is iterated to self-consistently determine both the smoothing 
length and the density according to the relation 

/ \ 1/2 

ft = (^^ j , (23) 

where rj is a constant relating the smoothing length to the 
average particle spacing. The procedure for doing this is 
described in IPrice fc Monaehanl tOOj] and IPricel J2004) . 
Enforcing this relation tightly has been found to signifi- 
cantly sharpen the resolution of a typical SPH simulation 
and in conjunction with the additional terms in the momen- 
tum (and energy) equation(s), in general leads to substan- 
tial improvements in accuracy (although perhaps, as we will 
demonstrate in this paper, with a loss of robustness). 

The pressure is calculated directly from the density us- 
ing a polytropic equ ation of state P = Kp^ and the standard 
cubic spline kernel llMonaghanlll99 is used. 

Damping towards the equilibrium solutions is achieved 
by applying a fo rm of the SPH a rtificial viscosity used for 
shock capturing llMonaghanlll997r) together with a damping 
in the force equation which is independent of resolution, 
given by 

^ = -0.02V + f, (24) 

where f is the force per unit mass. Note that since the pres- 
sure is calculated directly from the density, the kinetic en- 
ergy removed by the artificial viscosity and damping terms 
is not deposited as thermal energy but rather removed from 
the system. 

In the linear and non-linear oscillation solutions to 
be described, a small amount of dissipation is applied se- 
lectively using the artificia l viscosity switch described by 
iMorris fc MonaghanI il997l) . This switch turns on the vis- 
cosity terms in response to the magnitude of any negative 
divergence (ie. convergence) in the velocity field. In the os- 
cillation solutions this artificial viscosity is applied only to 
approaching particles. 

The time integration is achieved using a m anifestly re- 
versib le, second order integrator described by iMonaghanI 
i2005l) which in particular is reversible in the case of a vari- 
able step size (as for example when the step size is deter- 
mined from the Courant condition). At the present time, 
the reversibility condition only holds with the density cal- 
culated by direct summation and the pressure determined 
directly from the density, which is sufficient for the calcula- 
tions presented here (although it is fairly straightforward to 
generalise the integrator to the case in which the continuity 
and energy equations are also evolved). The differences be- 
tween using the reversible integrator and a simple predictor 
corrector method as commonly used in SPH are found to 
be minor, although in the course of the Toy Star tests we 
have found several aspects of the reversible integrator which 
must be treated with caution. 



3.2 Static structure 

The simplest test case for the two dimensional Toy Star is 
to verify the static structure. This can be done in one of two 
ways: using equal mass particles which are damped into an 
equilibrium configuration or by varying the particle masses 
according to the equilibrium density profile. In order to test 




Figure 1. Toy Star static structure. We place 1000 equal mass 
SPH particles on a lattice of square cells within the circle of unit 
radius and allow them to evolve under the influence of the linear 
force. The SPH particles are shown by the solid points after damp- 
ing to an equilibrium distribution whilst the solid circle shows 
the radius of the exact solution. The particles adopt a hexagonal 
lattice arrangement in the central regions and a ring-like arrange- 
ment at the outer edges 




Figure 2. Density profile in the Toy Star static structure with 
equal mass particles. The SPH particles are shown by the solid 
points after damping to an equilibrium distribution. The exact 
quadratic (p = 1 — r^) solution is given by the solid line. 



the robustness of our algorithm we examine both methods. 
The calculations are performed using an equation of state 



3.2.1 Equal mass particles 

In the equal mass particle case, we construct the initial con- 
ditions by placing the particles on a lattice of square cells 
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which is cropped to retain only those particles with a radius 
less than unity (chosen since it is the equilibrium radius of 
the Toy Star solution). Using a lattice spacing of 0.05, this 
results in 1256 particles. The particles are perturbed from 
the lattice with a random amplitude of up to | of the initial 
lattice spacing to remove any residual effects from the initial 
regular particle arrangement in the equilibrium solution. 

The system is then allowed to collapse under the influ- 
ence of the (axisymmetric) Toy Star force and damped until 
an equilibrium is reached (where typically we damp until 
Ekin ~ W~^° Etot)- The equilibrium particle configuration 
is shown in Figure Q (where for comparison the solid circle 
on this plot marks the exact solution radius). The particles 
can be seen to adopt a reasonably isotropic close-packed 
hexagonal lattice arrangement in the central regions of the 
Toy Star, whilst at the edge tend to form discrete rings. 
Experimenting with different initial particle configurations 
(including a completely random configuration) and different 
particle numbers result in the same distinguishing features 
(although with a particular caveat which is discussed below) . 
The configuration adopted by the particles is in general re- 
lated to the shape of the interpolation kernel used in the 
calculations. Thus the isotropic nature of the particle ar- 
rangement in the central regions results from the isotropy 
of the cubic spline kernel. The density profile of the equi- 
librium configuration is shown in Figure |5| and agrees well 
with the exact solution given by the solid line. The discrete 
rings of particles formed at the outer edges are also visible. 

The equal particle mass also provides a useful illus- 
tration of the 'particle pairing' instability which occurs for 
larger smoothing lengths, specifically where the constant of 
proportionality in 12311 . rj > 1.2 (where rj > 1.5 is found to 
be particularly unstable). In these cases the particles tend 
to initially form pairs, which then begin to coalesce. Even- 
tually these pairs clump together at the same location, with 
the resulting lattice adopting a close packed configuration 
similar to that shown in Figure (albeit at half the reso- 
lution since two particles effectively become one). This is a 
well-known defect of the cubic spline kernel, resulting from 
the fa ct that the force tends to zero at the origin (see f or ex- 
ample iThomas fc CouchmanllTooS ISwegle et aDll995ft . Us- 
ing ?7 > 1.5 results in the closest neighbour being placed on 
the turning point of the force curve (in the isotropic, equi- 
librium configuration), thus moving inwards (and therefore 
clumping) as the particles are compressed. 

We have also inves tigated this instability using the stan- 
dard SPH formalisms jHernauist fc Kcit3ll989l iBenz et alJ 
Il990l) whereby a simple average of either the particle 
smoothing lengths or the kernel gradients is used in the 
equation of motion. In these cases particles are still observed 
to clump, although the lattice configuration tends to retain 
defects due to the averaging procedure. With the new vari- 
able smoothing length formalism the particles thus appear 
to settle more readily into minimum energy states because 
of the smoother lattice configuration (Figure 0. This means 
that the pairing instability is slightly more pronounced in 
this case (although we have not made a quantitative assess- 
ment). 

A solution to this instability may be easily obtained by 
using a kernel where the force does not decrease towards the 
origin. However, the density estimate using such k ernels is 
in general quite poor. iThomas fc CouchmanI il992l) suggest 



a compromise approach in which the cubic spline kernel gra- 
dient is modified to remain constant in the region between 
the usual minimum and the origin whilst the unmodified 
cubic spline kernel is used in the density evaluation. Whilst 
this approach seems promising, it also introduces additional 
problems (for example, energy is not conserved exactly and 
the kernel gradient is no longer normalised). In this paper 
we simply use rj = 1.2 in all of the simulations, thus avoid- 
ing the instability, although it is our intention to investigate 
the issue further elsewhere (in particular to examine alter- 
native kernels to the cubic spline which have finite gradients 
at the origin but which do not lead to a significant decrease 
in interpolation accuracy). 

3.2.2 Unequal particle masses 

In the second case, the particle masses were varied accord- 
ing to the equilibrium density profile. As a first attempt 
the particles were arranged on a lattice of square cells, with 
masses corresponding to the density profile at their radial 
position. To avoid problems with very low mass particles at 
the outer edge we excluded particles within half of the initial 
lattice spacing of the equilibrium radius (r = 1 in this case) . 
Allowing the particle distribution to evolve under the infiu- 
ence of the linear force, with damping applied as described 
above, the particles at the outer edges of the configuration 
at first shift slightly in order to give a more circular edge 
(rather than the stepped edge given by the cropped lattice). 
When using the iterated smoothing length update the par- 
ticles were then observed to gradually shift from the square 
lattice and attempt to adopt a hexagonal- type arrangement. 
An equilibrium configuration which is mutually satisfying 
for all of the particles then becomes very difficult due to the 
unequal particle masses. The kinetic energy of the system 
is initially damped (to around W~^Etot) but climbs at later 
times (to around 10~^Etot) as the particle disorder spreads 
through the system. The particles at this stage appear to 
continually exchange between shifting blocks which do not 
seem able to settle to a minimum energy state (appearing 
not unlike the movement of shifting ice-fioes). Our under- 
standing of this phenomena is that the square lattice con- 
figuration with r] = 1.2 in II23I I is a quasi-stable state for the 
particles. This can be show n by a straightforward stability 
analysis in two dimensions (Morris "l99ef: 'B orve et al.ll2004l) 
which shows that the square lattice is unstable at certain 
ratios of the smoothing length to the particle spacing (al- 
though stable at = 1.2). Given enough sensitivity to their 
configuration (through the smoothing length iterations) the 
particles will shift into a more isotropic arrangement similar 
to that observed in the equal mass case, however in this case 
a minimum energy configuration is more difficult to obtain 
with unequal particle masses, since the masses have been 
set to correspond to the equilibrium density profile using 
the initial particle configuration. 

These results show that it is preferable to place the 
unequal mass particles initially on a lattice of hexagonal 
cells since this is the configuration to which they tend to 
evolve. 

We therefore set up the particles in this manner, with 
the (hexagonal) lattice centred on the origin and again with 
the particle masses set according to the equilibrium density 
profile. 1045 particles are used. In this case the lattice is able 
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Figure 3. Toy Star static structure with unequal particle masses. 
1045 SPH particles are set up in a hexagonal close packed lattice 
arrangement with their masses set according to the equilibrium 
Toy Star density profile. The particles are then allowed evolve un- 
der the influence of the linear force. The SPH particles are shown 
by the solid points after damping to an equilibrium distribution, 
whilst the solid circle shows the radius of the exact solution. 




Figure 4. Density profile in the Toy Star static structure with 
1045 unequal mass particles. The SPH particles are shown by the 
solid points after damping to an equilibrium distribution. The 
exact quadratic {p = 1 — r^) solution is given by the solid line. 



to damp to an equilibrium state since the particles do not 
shift from the lattice apart from a small initial adjustment 
at the outer edges (note that we also use a pmin in 12311 
as will be discussed in H^.^.l^ . The equilibrium distribution 
in this case is shown in Figure |3 and the density profile is 
shown in Figure |1| 



3.2.3 Reversible time integration 

The relaxation procedure for the Toy Star also revealed some 
pertinent aspects of the reversible integration algorithm. In 
particular if the timestep changes rapidly between steps, the 
energy in the reversible integrator could exhibit large oscilla- 
tions and under some circumstances become unstable. This 
is a result of the averaging procedure used to determine the 
timestep. In this paper we have used the arithmetic average 

5t^^^ = ^(5t° + 5t^), (25) 

where the superscripts 0, ^ and 1 refer to the stepsize calcu- 
lated at the beginning, middle and end of a given timestep. 
This expression is rearranged in order to determine the 
next timestep St^ from the previous timestep 5t° and the 
timestep determined from the Courant condition at the half 
step St^^^. The pitfall of this method is that if the timestep 
changes rapidly (for example, from a very short step to a 
very long step) then the timestep begins to oscillate between 
two values (short-long-short-long). Also there is nothing to 
prevent the timestep 5t^ from taking a negative value. In 
this situation the time advance is achieved by taking a long 
step forward and a slightly shorter step back to achieve a 
small net forward step. Whilst the specific problem of neg- 
ative timesteps can be fixed by using a different averaging 
procedure (for example using the geometric mean), we have 
found that the oscillations can also cause instability if St^ 
greatly exceeds the Courant condition. This can occur be- 
cause a very short St*^ implies a very long 5t^ (where by def- 
inition St^''^ meets the Courant condition). If the Courant 
constraint also changes between steps the result can be that 
particles are evolved on a timestep which is too large to cap- 
ture the change in physical quantities accurately, resulting 
in instability. 

In practice these oscillations can be largely avoided pro- 
vided that some care is taken in the setup and evolution. 
For example in the initial runs we began the integration 
as we had done for the predictor-corrector method using a 
timestep of zero (which simply runs through the force evalu- 
ation without evolving the particles). In the reversible inte- 
grator this immediately leads to timestep oscillations since 
the zero timestep is then used in the averaging procedure. A 
second example is that it is common practice to reduce the 
timestep before output dumps in order to reach the spec- 
ified time of output exactly. This means taking a shorter 
step just before an output dump but for the reversible inte- 
grator this will again trigger the timestep oscillations. Thus 
an interpolation (rather than an evolution) must be used 
to output the data at a specific time. Finally, in order to 
prevent any potential instability we place a check in the 
timestepping algorithm which resynchronises St'^, St^^^ and 
St^ to the current value of St^^^ if St^ exceeds the Courant 
condition (but without a safety factor). This means using a 
reasonably low safety factor for the Courant condition used 
to find St^'''^, in order to allow sufficient variation between 
St'\ 5t^^^ and 5t^. Obviously any timestep resychronisation 
means that the evolution is no longer globally reversible, 
although the sections of the evolution between resynchroni- 
sations will be reversible. 

It should be noted that methods do now exist for re- 
versible integration of SPH-type systems which do not suf- 
fer from timestep oscillation instabilities. These methods in- 
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Figure 5. Axisymmetric (s=0) modes j=2, 4, 6, 8, 10, 12, 14, 16 and 18 after 1 (theoretical) oscillation period using 1045 particles with 
masses initially varied to give the equilibrium density profile. The axisymmetric modes exhibit strong oscillations near the outer edges 
of the Toy Star which can be difficult to capture numerically. 




Figure 6. Particle distribution in the unequal particle mass run 
after 10 oscillations of the axisymmetric j = 2 mode. In the un- 
equal mass case using h = ri{m/ p)^/'^ can lead to large smoothing 
length gradients at the edge of the Toy Star, causing some disrup- 
tion at the outer edge (left figure). Using h = ri[m/(p + pmin)]^^^ 
(right figure) prevents this problem from occurring. 



volve a rescaling of the time variable so that the timestep is 
a continuously variable quantity (Leimkuhler, private com- 
munication). The application of such methods to SPH will 
be investigated elsewhere. 



3.3 Linear Oscillations 

The linear oscillations of the two dimensional Toy Star may 
be examined by applying a velocity perturbation to the equi- 
librium configuration. We examine both the equal mass and 
unequal particle mass cases for both the axisymmetric and 
non-axisymmetric modes. The perturbation in each case is 
given to the velocity components according to I22II (where 
the higher order modes are calculated using the recursion 
relation ilUt ) and normalised such that the total kinetic en- 
ergy 



i(0.05c,o)', 



(26) 



where Cso is the sound speed at the origin. 

The SPH simulations in each case were run for 10 peri- 
ods of the given mode. The frequencies of oscillation in the 
numerical solutions may then be obtained from the kinetic 
energy evolution. The simplest method is to measure the 
spacing of minima and maxima in the kinetic energy. We 
use this method as well as calculating the full power spec- 
trum of the kinetic energy evolution using the Lomb/Scargle 
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Figure 7. Density perturbations in the axisymmetric (s=0) modes j=2,4,6,8,10,12,14,16 and 18 using 10053 unequal mass particles. The 
difference between the density and the equilibrium density profile is plotted for each mode after 1.25 oscillation periods. 
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Figure 8. Density perturbations in the j=0 modes with 8=2,4,6,8,10,12,11,16 and 18 using 10053 unequal mass particles. The difference 
between the density and the equilibrium density profile is plotted for each mode after 1.25 oscillation periods. 
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Figure 9. Time evolution of tlie first few axisymmetric (s=0) modes (from top to bottom j=2,4,6,8 and 10) using 1000 equal mass 
particles. The left hand side shows the time evolution of the kinetic energy over 10 periods for each mode, whilst the right hand side 
shows the dominant frequencies of oscillation from the power spectrum (periodogram). Prom this plot we can see that the modes up to 
j = 8 are captured well at this resolution (although the amplitudes are somewhat damped), whilst the higher order modes (in this case 
j > 10) tend to decay into their lower order counterparts. Note that two oscillations of the kinetic energy correspond to one oscillation 
period. 



periodoRr am for un evenly sampled data jLombll976l:IScargld 
Il982) (see iPress et al...l99a) . 



3.3.1 Axisymmetric modes 

The linear, axisymmetric modes of the Toy Star are shown 
in Figurel^l where the numerical solution in the unequal par- 
ticle mass case is plotted for each mode after 1 oscillation 
period. The perturbation solution in each case is given by 
the solid line. The modes change rapidly near the boundary 
because the low density material respond to the push by 
getting a large velocity. This shows up in the rapid change 
of the hypergeometric functions which describe this varia- 
tion. However, these strong oscillations can be difficult to 
capture numerically. In the equal mass particle case at this 
resolution the modes are significantly damped by the artifi- 
cial viscosity terms which act preferentially on the outer re- 



gions of the Toy Star^. In the unequal particle mass case the 
large fluctuations at the outer edges are better resolved since 
the resolution is constant throughout the Toy Star, however 
some problems occur at the edges due to mixing between 
particles of different masses caused by a combination of the 
large smoothing length gradient at the outer edge, the ve- 
locity ffuctuations from the oscillations and the presence of 
very low mass particles there (Figure [S] left panel). 

The large smoothing length gradient at the edge of the 
Toy Star is largely a defect of our use of 1231 in the case 
of unequal mass particles. Whilst the problem also occurs 
using standard SPH formalisms, it is accentuated by our 
use of the variable smoothing length formalism since the 



^ Note that the SPH artificial viscosity has a kinematic viscosity 
coefficient which has the form oc (hCs) which is constant near 
the edge since h oc l/-y/p and Cs oc y/p (when P oc p^). Whilst 
for other adiabatic equations of state the coefficient becomes infi- 
nite, in this case the damping is large in the outer regions purely 
because the second derivatives of the velocity are very large. 



large gradients in the smoothing length are explicitly incor- 
porated into the force terms. The relation I23II may cause 
problems in this case because, ideally, the smoothing length 
should relate to the particle number density (or roughly, the 
number of neighbours) rather than the mass density. In the 
case of equal mass particles ()23|l indeed has this effect. The 
problem in the unequal mass Toy Star is that where the 
smoothing length should increase at the outer radius (due 
to the decrease in particle number density), in practise the 
smoothing length in fact decreases at the edge due to the 
drop-off in particle mass affecting the numerator in 1231 . A 
simple fix to this problem is to modify the relation using a 
density floor, giving 

/ \ 1/2 

V P + Pmin / 

where pmin = min(ppreu), where ppmv is the density on the 
particles from the previous timestep. This modiflcation pre- 
vents the smoothing length gradients from becoming too 
large at the edge. 

Using a density floor is not particularly desirable in 
general. In fact it is quite simple to derive a gener- 
alisation of the variable s moo thing length form alism of 
ISprineel fc Hernauisd ll2002h and lMonaghanl (jiool) in which 
the smoothing length is a function of the particle number 
density rather than the mass density^. However, since this 
is a paper on Toy Stars rather than on SPH formalisms the 
detailed description of this formalism will be presented else- 
where and we simply use (1271 in this paper. In terms of the 
frequency estimate there is very little difference between the 
results using H23|l and those using the density floor. In fact 
the disruption seen in ® is largely cosmetic for short simu- 
lations since the very low mass particles do not have a strong 
influence on the overall evolution. For longer simulations the 
effect is more problematic as the disruption in the particle 
distribution spreads to the inner regions. It should also be 
noted however that these effects also become less significant 
as the total number of particles is increased. 

The kinetic energy evolution in the first 5 axisymmet- 
ric modes in the case of equal mass particles is shown in 
Figure El The left panel shows the kinetic energy evolution 
whilst the right panel shows the power spectrum of the left 
panel, showing the dominant frequency in each case. The 
modes show significant damping although the dominant fre- 
quency is captured well by the numerical solution up to the 
j — 10 mode which can be seen to decay rapidly into the 
j — 8 mode. 

The frequencies in the numerical solutions are shown in 
Figure ITUl The top panel shows the absolute frequency for 
each axisymmetric mode plotted as a function of the mode 
number j. The frequency plotted in each case is the dom- 
inant frequency measured from the power spectrum of the 
kinetic energy evolution for the simulation of that partic- 
ular mode. The exact frequencies are given by the crosses, 
where a solid connecting line has been drawn for clarity. 
Results using equal mass particles are given by the open 

^ A similar formalism was desc ribed bvlO tt fc SchnetteJ i200,'j) 
and derived self-consistently in iPricj but lias the disad- 

vantage of also modifying the density evaluation which we have 
found to be somewhat problematic. 
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Figure 10. Frequencies of the axisymmetric modes at a reso- 
lution of 1000 particles. The top panel shows the absolute fre- 



quencies whilst the bottom panel shows the fractional error in 
the numerical solution. The simulations using unequal particle 
masses, denoted by filled circles (a connecting line (dashed) has 
been plotted for clarity) are within 2 percent of the correct fre- 
quencies for modes up to j = 20. The results using equal mass 
particles (open circles, dotted line) show errors of up to 10% for 
modes above j = 6 although the general trend is still observed. 




^ -0,1 - 



5 10 15 20 

phi mode (s) 

Figure 11. Frequencies of the j = modes at a resolution of 1000 
particles. The top panel shows the absolute frequencies whilst the 
bottom panel shows the fractional error in the numerical solution. 
The simulations using unequal particle masses, denoted by filled 
circles (where again a connecting line (dashed) has been plot- 
ted for clarity) are within 2 percent of the correct frequencies 
for modes up to s = 16. The results using equal mass particles 
(open circles, dotted line) show somewhat lower errors but only 
for modes up to s = 8 above which the modes are largely damped 
out. 



circles (with a dashed connecting line for clarity) whilst the 
unequal particle mass results are indicated by the filled cir- 
cles (with a dotted connecting line). The lower panel shows 
the error in the numerical frequency (shown as a fractional 
deviation from the exact value) in each case. The unequal 
particle mass simulations at this resolution show errors of 
< 2% for all modes up to j = 18. The results using equal 
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Figure 12. Density perturbations for tiie mixed modes using 10053 unequal mass particles. From left to right the radial mode number 
varies from j = 2 (leftmost column) to j = 12 (rightmost column) whilst from top to bottom the angular mode number varies from s = 2 
(top row) to s = 12 (bottom row). The difference between the density and the equilibrium density profile is plotted for each mode after 
1.25 oscillation periods. 



mass particles have comparable accuracy up to the j = 6 
mode but thereafter show errors of ~ 10%. One interest- 
ing point to note in the equal particle mass case is that the 
j = 10 mode seems to decay to the j ~ 8 mode and corre- 
spondingly the j = 20 mode seems to decay to the j = 16 
mode. The frequencies obtained with the time-reversible in- 
tegrator are indistinguishable from those obtained using the 
predictor-corrector method, provided that the timestep does 
not change rapidly between timesteps as discussed in H3.2I 



3.4 Non-axisymmetric oscillations 

3.4-1 j ~ modes 

The two dimensional Toy Star can also oscillate non- 
axisymmetrically. The frequencies in the numerical solutions 
at a resolution of 1000 particles for the j = modes (j=0, 
s=2,4,6 etc) are shown in Figure 1111 As in Figure 1101 the 
absolute frequencies are shown in the top panel whilst the 
fractional error is shown in the bottom panel. Using unequal 
mass particles (filled circles, dashed line) the frequencies in 
the numerical solutions show good agreement with the ex- 
act solutions (errors < 2%) for modes up to s = 16, above 
which the modes are damped out. The results using equal 
mass particles (open circles, dotted line) show somewhat 
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lower errors (< 0.6%) but only the modes up to s = 8 are 
captured in this case. 

3.4.2 Mixed modes 

Finally the Toy Star can oscillate with modes which have 
both J ^ and s 7^ 0. The density perturbation in the 
modes with 2 < j < 12 and 2 < s < 12 is shown in Figure^] 
The numerical frequencies of the perturbations show similar 
accuracy to the symmetric modes. 



4 EXACT, NON-LINEAR SOLUTIONS 

As mentioned earlier, one of the most useful features of the 
Toy Star is that it is possible to get exact, or very accurate 
solutions of the full non linear equations. These provide ex- 
cellent tests of the ability of a numerical code to follow a 
surface into a region devoid of matter. We assume that the 
velocity field is a linear function of the cartesian coordinates 



,,11 , ,,12 
Vx = V x + V y, 

Vy = V'^^x + V^^y, 



(28) 
(29) 



and that the density is given by the following second degree 
function of the coordinates 



H{t) - [C{t)x'^ + 2B{t)xy + D{t)y\ 



4.1 Axisymmetric nonlinear solutions 



For this case we choose V 

r21 



a{t) and 



-V'^' = ~f3{t) in which case the velocity field can be written 

V = Q(t)r + /3(t)z X r, (31) 

where the first term produces an axisymmetric expansion 
or contraction, and the second term gives a rigid rotation 
with angular velocity /3(t)z where z is a unit vector in the 
z direction. The motion takes place so that, at any time, 
each element of fluid has the same angular velocity. The 
reader will appreciate that the solution has many features 
in common with the symmetric pulsation of a polytropic 
star, for example a model White Dwarf, under gravity. 
For this axisymmetric case C — D, B = 0, so that 



9^-' = H{t) - c{ty. 



(32) 



Substitution of the assumed forms for v and p into the 
acceleration and continuity equations, with P = Kp'' and 
replacing the Lagrangian derivative by the Eulerian accord- 
ing to 

and equating powers of x and y we get 



a = -a + f3 -\ ; SZ , (34) 

7-1 

$ = --2a/3, (35) 
H = -2a(7-l)7i", (36) 
C = -2Ca7. (37) 
Where the a denotes the time derivative of any function a. 



If 7 = 2, a further derivative of 13411 with respect to 
time using 13511 to remove the derivative of /3 results in the 
following equation for a 



0. 



(38) 



If a is sufficiently small this equation can be approximated 
by 



a + AaQ, = 0. 



(39) 



with solution a = Asin{2^lt + c) where c is an arbitrary 
constant. An exact solution of the full non linear equation 
can also be found. The analysis is given in Appendix IHl 

The set of ordinary differential equations 13411 - 137|l can 
be integrated with high accuracy using any standard stable 
ODE integrator. 



4.1.1 Conserved quantities 
The total mass is 

■kHT^ (7 - 1) 



(30) 'HL 



and must be conserved by the set of equations 1341 - 137(1 . To 
see that this is the case we eliminate a from the last two of 
these equations to get 

dH _ {'y-l)H 
dC ^ ^ ' 

from which H'^ oc C"^ showing that A/ is constant. The 
angular momentum Jz should also be conserved by the equa- 
tions. Jz is given by 



(41) 



J, = 2-KI3 



p{xvx-yvy)dxdy = 2tvI3 / [//-Cr^]~r^dr. (42) 



(43) 



Jz = 



Completing the integration shows that 



C2 27(27 - 1) • 

But from 13511 and 13711 we find oc C. Substitution into 
H43II and noting that oc C^~^ we find Jz is constant. Of 
course the system must conserve AI and J^ and the foregoing 
analysis merely confirms that the equations we have derived 
are correct. 



4.2 The Lagrangian description 

In this section we consider the solutions in terms of the 
change in the functions H, C and /3 for a given element of 
fiuid. Defining the function S{t) by 



S{t) = exp{- I adt). 
Jo 



we can write 



p = 


f3oS\ 


H = 




C = 





and from the continuity equation © 
P = poo 



(44) 

(45) 
(46) 
(47) 

(48) 
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where we assume the motion begins at t = 0. These equa- 
tions can be interpreted in the following Lagrangian way. 
For any element of fluid (particle), the initial value of one 
of these functions, say P, is denoted by /3o- The equation 
for j3 in terms of S gives its value for that particle at any 
time. Similarly for the other quantities. This also enables 
one quantity to be written in terms of the others. For exam- 
ple we deduce H(5 oc C. The x, y coordinates of a particle 
can be obtained from the velocity according to 

J = c^^-Py, (49) 

dy 
dt 



ay + I3x, 



(50) 



From these equations we easily deduce that the radial coor- 
dinate r and angular coordinate of the particle are given 
by 



I3{t)dt, 



(51) 
(52) 



from which we deduce that for any particle Cr^ cc H as \a 
required by the density equation (I32II and 14811 . 

A further useful result is that equations 1341 and 13511 
can be combined together with the relation C = to 
give 



da 
dp 



2a/3 



a^ + 0.'^ 



13' 



2'yK 
7-1 



a/3' 



(53) 



This equation can be integrated exactly taking a as a new 
dependent variable to give a linear differential equation. In 
this way we find 

2it"7C 



(7-1)^ 



+ kf3, 



(54) 



where k is an integration constant. This equation can also be 
derived from the constancy of energy. 15411 defines a closed 
curve in the a (3 plane showing that the motion is a periodic 
oscillation (demonstrated in Figure ITSl . 

4.3 SPH simulations of the non linear 
axisymmetric motion 

The non-linear axisymmetric Toy Star is set up by perturb- 
ing the equilibrium solution with a velocity of the form 13111 
where the parameters a(0) and /3(0) specify the amplitudes 
of the initial expansion and rotation respectively. 

Since the axisymmetric mode given by 13111 is a non- 
linear, exact solution the choice of amplitudes is completely 
arbitrary. We chose for simplicity, ao = 1 and = it. This 
means that the expansion velocity is everywhere supersonic 
(ie. highly non-linear) and that initially the rotation rate 
was equal to the oscillation period (note that the actual ro- 
tation period during the evolution depends on the periodic 
manner in which (3 varies). The SPH simulation was run 
by perturbing the equal mass particle equilibrium configu- 
ration. No artificial viscosity was applied during the evolu- 
tion. The density and radial velocity are shown in Figure [T^ 
after 100 oscillation periods. The agreement with the exact 
solution after 100 periods is extremely good, with both the 
density and velocity profiles maintained almost exactly. 

The time evolution of the total kinetic energy in the 




Figure 13. Density and radial velocity in the non-linear Toy 
Star solution after 100 oscillation periods. The SPH particles are 
indicated by the solid points whilst the exact solution is given by 
the solid line. 




Figure 14. Time evolution of the nonlinear axisymmetric mode 
of the Toy Star using 1000 equal mass particles, where initially 
= 1 and /3 = TT. For clarity only the first 20 oscillation periods 
are plotted. The amplitude is maintained exactly by the SPH 
solution for more than 1000 periods. 



simulation is shown in Figure ITU showing only the first 20 
periods for clarity. The amplitude of the oscillation is main- 
tained exactly by the SPH solution. 

Finally it is useful to plot the evolution of the La- 
grangian quantities a and /? in the numerical solution. In 
order to do so we calculate Qi = (v • f)i and /3i = (f x 'v)^^i 
for each particle i after every timestep and then compute 
the average over all the particles. These average values of a 
and 13 are plotted in Figure [131 where each point in the plot 
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Figure 15. The Lagrangian quantities a and /3 measured from 
tlie SPH solution using 1000 equal mass particles. The average 
values of a and /3 on the particles has been plotted every timestep 
for over 1000 oscillation periods {t=3200, 72,863 timesteps). 



corresponds to the values at a given timestep. The curve 
agrees exactly with the analytic solution given by (I54II . 

The results of this simulation demonstrate that SPH has 
extremely good conservation properties and can maintain 
an accurate evolution of a non-linear system over long time 
intervals. We have experimented with a range of values for 
Qo and Po all of which show similar results. In practise the 
necessary addition of artificial viscosity into SPH algorithms 
in order to handle shocks leads to unphysical damping of 
self-similar motion such as that which occurs in the non- 
linear Toy Star. This test would therefore be an excellent 
benchmark for switches designed to turn off the artificial 
viscosity away from shocks where it is not needed. 



4.4 Non-axisymmetric solutions 

If the expressions H28|l - H3()| l for the velocity and are 
substituted into the equations of motion and the coefficients 
of powers of x and y and xy the following equations are 
obtained from the acceleration equation 





2'yK 


dt 


7-1 


dl/22 


27A' 


dt 


7-1 


dV'^ 


27A' 


dt 


7-1 


dV^i 


2'yK 


dt 


7-1 



c - (v")^ - y^V^^ 



D - (v^^f - v^v^^ - n^ 



and from the continuity equation we get 



(55) 
(56) 
(57) 
(58) 

(59) 



^ = -2Cy"-(7-l)C(V''+\/^^)-2BV'\ (60) 
at 

^ = -2W^'' - (7- 1)D(V" +y'') -2BV'^ (61) 
dt 



dB_ 

dt 



These 8 ordinary differential equations describe the system. 
They can be integrated as accurately as desired using stan- 
dard methods for differential equations. 

4.4.1 Conserved Quantities 
The mass is given by 

M=l I {H - {Cx^ + 2Bxy + Dy'^))^^'^"'-'^Uxdy (63) 



This integral can be evaluated by rotating the coordinate 
system so that the quadratic form in the integrand is reduced 
to a sum of squares. The final result is 



M = n 



7 - 



l\ 



7 J VCD - 52 
From the differential equations we find 
d{CD - 



(64) 



dt 



= -27(CL» - B^)(V'" + V^^), 



(65) 



which, combined with 1)59^ shows that A/ is conserved. 
The angular momentum is given by 



Jz = 



{H - {Cx^ + 2Bxy + Dy'^)f'^''-^'>Xdxdy, (66) 
where we have defined 

X = (x V + xy{V^^ - V") - y V') . (67) 

Rotating the coordinate system as before the integral re- 
duces to 



_ ^7-1)^ H(2^-l)/(7-l) 
7(27 - 1) (CD - £2)3/2 ' 

where we define 

T = V^^D - V^^C - B{V'^^ - y"). 



(68) 



(69) 



From the differential equations above it is possible to show 
that the quantity T satisfies the equation 



^ = -(7+l)T(V" 4- V"") 

In terms of the variable T defined by 



T = exp[~ I {V^^ + V^^)dt 



(70) 



(71) 



(62) 



we can write from 159L (1651 and the previous equation 

H = HqT'-''-^\ (72) 
(CD-B^) = (CD -5^)07^", (73) 
T = ToT^+V (74) 

When these results are substituted into the previous expres- 
sion for Jz it is found to be constant, as expected. As in the 
axisymmetric case proving the constancy of AI and Jz from 
our set of differential equations merely confirms that they 
are correct. 

Another conserved quantity is the circulation. The cir- 
culation around any loop for this two dimensional fluid is 
given by 

C = ^v-dl = J j V X vdxdy. (75) 
This is a constant of the motion so that its time derivative 
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Figure 16. Particle distribution during the evolution of the nonlinear, non-axisymmetric mode of the Toy Star using 1000 equal mass 
particles, where initially Vii = 1, V12 = 1/2, V21 = 1 and V22 = 1/4. A small damping of the amplitude may be observed in the SPH 
solution due to the artificial viscosity applied to prevent particle interpenetration during the compression phase. 



must vanish. Noting that V x v = {V'^^ — V^^)z and that 
from the equations of motion 

_ yl2 ^ ^y21 _ y-12)yy^ (77) 

and 

P = poT, (78) 

we see that (V x v)/p is unchanged following the motion 
and C is conserved. 

4.5 SPH simulations of the non linear non 
axisymmetric motion 

The non-axisymmetric exact solutions present a more de- 
manding test problem since the Toy Star changes shape 
throughout the oscillation. This changing shape must then 
be followed through the expansion and contraction phases, 



maintaining the free boundary. The non-axisymmetric Toy 
Star therefore presents an ideal problem for numerical codes 
designed to simulate astrophysical systems with changing 
geometries and free boundaries. 

The simulation is set up by again perturbing the equi- 
librium solution, in this case with a velocity of the form 
Il28l l- ll29| l The four parameters Vii,Vi2,V2i and V22 specify 
the amplitude and geometry of the perturbation. Due to the 
changing shape the SPH solution in this case requires some 
artificial viscosity in order to prevent particle interpenetra- 
tion during the compression phase. We apply this using the 
iMorris fc Monaehanl lll997fl switch as discussed in H3.ll 

The particle distribution during the evolution of the ex- 
act, non-axisymmetric mode is shown in Figure [TSl where 
initially Vii = 1, V12 = 1/2, V21 = 1 and V22 = 1/4. For 
comparison with the exact solution we solve Il55t - 162j us- 
ing a simple second order modified Euler predictor-corrector 
method, using the conservation of mass (1641 as a check on 
the quality of the integration. The solid line shown in Fig- 
ure ITSlis the curve corresponding to the edge of the Toy 
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Star (ie. where p = 0) at the appropriate times. The SPH 
particles adjust to the changing shape quite well apart from 
a damping of the amplitude with time caused by the appli- 
cation of artificial viscosity. We have performed a range of 
simulations using different values of the initial parameters 
which in general show very similar results. 

For simulations with very strong compression the parti- 
cles can clump together in the manner described in H3.2l due 
to the force being zero at the origin of the cubic spline kernel 
used in the calculations, despite using h — 1.2(m/p)^/^. This 
is because a strong compression can push the particles neigh- 
bours close enough to be in the region of the kernel where 
the force decreases towards the origin, causing a clumping 
instability for positive pressures in compression similar in its 
effect to the well known in stability for negative stresses in 
tension iSweele et al As discussed in >l3.2l the remedy 

for this is to use a kernel with non-zero derivative at the ori- 
gin, an investigation of which will be performed elsewhere. 



5 SUMMARY AND DISCUSSION 

The primary aim of this paper has been to provide a set 
of benchmarks for simulations of gaseous disks of the kind 
that arise in star formation. The systems, which we call Toy 
Stars, are similar to their astrophysical counterparts in that 
they consist of compressible gas, held together by an attrac- 
tive force which results in the gas having an outer surface 
where the density and pressure fall to zero. They therefore 
provide tests which are quite different to the usual tests 
based on flow in periodic rectangular regions. Furthermore, 
not only can the linear modes be found in terms of known 
functions, but non linear solutions can be found in terms of 
a small number of differential equations. 

The results described and discussed in this paper can 
be summarized as follows: 

(i) The static structure has been simulated using SPH 
with a variety of initial configurations using both equal mass 
and unequal mass particles. The calculations use approx- 
imately 1000 particles. With this number of particles the 
agreement of the static model with theory is very good. The 
final arrangement of the particles is different depending on 
whether they are equal mass or unequal mass. Recommen- 
dations concerning damping the particles to equilibrium are 
contained in the text. 

(ii) The SPH simulations use the equations which result 
from using a variational principles and taking the resolu- 
tion length to be a function of the density. The form of this 
relation, and effects that occur for different choices of pa- 
rameters is discussed in detail. 

(iii) The theoretical linear modes of oscillation have been 
derived and they have been calculated using SPH. Because 
the resolution of the SPH calculations is not high the sim- 
ulation would not be expected to recover modes high or- 
der modes. In fact, with 1000 particles the particle spacing 
is typically 0.06 (though larger near the edge in the case 
of equal mass particles). Accordingly, since several particle 
spacing between nodes is normally required to give reason- 
able accuracy, the linear modes of the Toy Star with mode 
number greater than about 3 would not be expected to be 
accurate. Remarkably the SPH simulations give good fre- 
quencies for radial and azimuthal mode numbers (s) much 



larger than this. The modes, however, are not simulated as 
accurately. We note, in particular, that in the outer layers of 
the Toy Stars where the density and speed of sound falls to 
zero, the modes vary rapidly, and they would only be repro- 
duced accurately by a simulation with many more particles 
(as for example in Figures 171 151 and 1121 used to illustrate the 
higher order modes). 

(iv) Our analysis of a class of non-linear modes reduces 
the partial differential equations to a set of 8 non linear or- 
dinary differential equations. These have been derived, with 
a different interpretation, in the context of the oscillation of 
water in basins. The equations describe a rotating oscillat- 
ing ellipse of gas which is similar to what would be expected 
when the attractive force is gravity. The SPH simulation of 
this system is in good agreement with the accurate integra- 
tion of the ordinary differential equations. 

(v) Our results show that SPH simulations of compact 
rotating masses of gas are very satisfactory. The actual ac- 
curacy depends, of course, on the resolution and therefore 
the number of particles used. 
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APPENDIX A: LINEAR MODES 
Al Axisymmetric modes 

For the axisymmetric modes s — 0. The lowest axisymmetric 
mode has j = 2 and setting 

we find 



(Al) 



7 



-X- 



For the next mode (j — 4) 
^ ^ \}~r'^ + 16(7-1) ^ 



V 



4cr 



2(7-1) rl 



(A2) 
(A3) 

(A4) 
(A5) 



where i/f — 4(47 ~ 2)/(7 — 1). Higher order terms can be 
obtained easily using the recurrence relation for the a^. 



A2 Non-axisymmetric modes 

The non axisymmetric case s = 1 is shift of the system. To 
see this in the simplest case we note that the lowest mode 
has j — and the spatial part of the perturbation rj for the 
choice sin 6* for the angular variation, has the form 



D — aor sin 9 = UoX, 



(A6) 
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so that the x component of the centre of mass becomes 



7-1 



77 ) xrdrdO. 



(A7) 



When s = 1, D is an odd function in x so the x compo- 
nent of the centre of mass is non zero. This mode therefore 
corresponds to the centre of mass shifting back and forth in 
the X direction which is impossible for an isolated object. 
On the other hand, it is entirely possible in a lake of water 
perturbed by tidal forces or a glass of wine moved back and 
forth. It can be shown that, this mode (s = 1 and j = G ) 
results in the surface of the wine, or the lake, remaining flat. 
If we choose the angular variation cos6 for the modes the 
shift in the centre of mass is in the y direction. All modes 
with s = 1 produce a shift in the centre of mass because 
they are odd in either x or y. For our astrophysical problem 
we therefore jettison the s = 1 modes as being unphysical. 

The first non axisymmetric mode which is physical has 
s = 2. The case j = has 



C(0 = ao - 



The next s = 2 mode has j = 2 and 



C(r) = aoX" (^1 



37-2 2 
3(7-1) 



(A8) 



(A9) 



with = 4(37 - l)/(7 - 1). For this case the velocity com- 
ponents are 



" a 3(7-1 r| 

V. - ^.fl-l^^ 
a \ 3(7-1 ri 



(AlO) 
(All) 



APPENDIX B: EXACT SOLUTION FOR NON 
LINEAR, AXISYMMETRIC EVOLUTION 

In this Appendix we give the outline of the solution of the 
equation 



Q + 6aa + 4a^ + 4Q.^a = 0. 



(Bl) 



which applies to the axi-symmetric, non-linear oscillation 
when 7 = 2. 

Let V — a and 



^a^. The equation becomes 



= -(40' + 6V + 8g). 

dq 



(B2) 



Define a new variable = 4f2' -|- 8q and substitute into the 
previous equation. The resulting homogeneous equation can 
be solved by setting V = FC, and solving for F. We find 



J2F + If ^ 
^ AF + l) 



(B3) 



where c is a constant. Returning to the variable V , and after 
some algebra, we get 



V = 



(B4) 



Recalling that = 4(n' + a^) we can write a in terms of 
^. Next we use a new variable 77 = y/c — in place of ^ and 
finally set x = ri + ^/c. This reduces the integral to the form 



dx 



(B5) 



which can be integrated and combined with the transformed 
IIB4I I to give 



x\/Ac 



16J72 



= sin {2D,{t + to). 



(B6) 



Transforming back to 77, solving the resulting equation for 
■q^, and then converting this to an equation for we get 



2 (o"C 
4q = ^ 



(B7) 



(B8) 



where cr = (c - 40.'^), and 6 = 2f7(t -I- to). Finally, 

_ Q. cos 9y^c - 4^2 

4- sin - 4^2 ' 

By a suitable choice of to the cosine and sine can be inter- 
changed. Furthermore, from IIB3II it can be shown that if a 
is small then \/c — 40? is small. In this case the solution 
reduces to that given after 
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